Turduckening black holes: an analytical and computational study 
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We provide a detailed analysis of several aspects of the turduckening technique for evolving 
black holes. At the analytical level we study the constraint propagation for a general family 
of BSSN-type formulation of Einstein's field equations and identify under what conditions the 
turducken procedure is rigorously justified and under what conditions constraint violations 
will propagate to the outside of the black holes. We present high-resolution spherically sym- 
metric studies which verify our analytical predictions. Then we present three-dimensional 
simulations of single distorted black holes using different variations of the turduckening 
method and also the puncture method. We study the effect that these different methods 
j—i have on the coordinate conditions, constraint violations, and extracted gravitational waves. 

I , We find that the waves agree up to small but non-vanishing differences, caused by escap- 

^jjy ing superluminal gauge modes. These differences become smaller with increasing detector 

location. 
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^ . I. INTRODUCTION 

On 

qq ' In a previous publication [1[ we discussed the turduckening approach to numerical simulations 

of black holes in Einstein's theory. The technique relies on initially smoothing the data inside 
each black hole and solving the Einstein evolution equations everywhere at later times. The idea 
was first proposed in Ref. 0] under the name of "free black hole evolution". It shares many 
similarities with the "stuffed black hole" 0, Q] and "magic matter" [^] approaches. In Ref. [l[ 
we presented a particular implementation that works in practice for binary black holes. We also 
provided justification for our implementation, and numerical evidence of the geometrical picture 
behind it. Complementary results based on very similar ideas were independently found and 
presented in Ref. 0] under the name "filling the holes" . 

The intuitive rationale behind the turduckening approach is that the physics in the exterior of 
a black hole should be causally disconnected from the unphysical smoothing in the interior. This 
is the same rationale behind black hole excision 0, 0|, but here one proceeds in a different way. 
In particular, one does not need to place an inner boundary per black hole in order to remove the 
interior. The computational domain in this technique is trivial from a topological point of view, 
and therefore the discretization remains simple. Thus, the method shares the simplicity of the 
moving punctures technique [jj, |l3j but is not restricted to puncture-type initial data and does not 
require regularization of the equations near special points. 
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In this paper we extend our analysis of the turducken technique, concentrating on both concep- 
tual and practical issues. 

We begin in Section [II] by describing the formulation of the equations that we used in Ref. [l|, 
which is a specific version of the BSSN-type family. We analyze in detail the hyperbolicity of 
both the main system and the subsidiary (constraint) system, placing particular emphasis on the 
propagation speeds of constraint violating modes. It is well known that in the Einstein equations 
the "true" degrees of freedom are coupled to coordinate and constrained degrees of freedom. One 
therefore needs to guarantee that, for the formulation of the Einstein evolution equations and the 
gauge conditions being used, the smoothing in the interior of each black hole does not affect the 
"physics" in the exterior. This is a non-trivial condition and, in fact, it is formulation and gauge 
dependent. In Section |TI] we show that there are some versions of the BSSN equations where this 
condition does not hold, and where constraint violations that originate in the interior of the black 
hole do propagate to the outside. However, we are also able to identify a class of BSSN-type 
equations for which we can rigorously guarantee that constraint violations inside the black hole do 
not leak to the outside. 

Next we concentrate on the issue of whether gauge modes can escape from the interior of the 
black hole. The gauge conditions that we use are those of the moving punctures technique. In 
Section |TI] we show that some of the characteristic speeds depend on the solution itself. Therefore it 
is not possible to determine a priori whether or not some modes will become superluminal. There 
is nothing wrong with modes leaking from the black hole interior, as long as these modes represent 
the gauge freedom inherent in the evolution problem. It is nevertheless of conceptual and practical 
importance to understand how the turduckening procedure might affect the gauge outside the black 
hole. Below we turn to this point by analyzing the numerical data. 

Having analyzed the system of equations at the continuum, and in particular, having shown that 
at that level the turduckening procedure does not introduce constraint violations to the exterior of 
a black hole, we address the discretization and numerical implementation in the following sections. 
We begin in Sec. IIIII with a brief description of the numerical codes that are used in this paper. 
In Sec. IIVI we evolve turduckened initial data for a Schwarzschild black hole with a spherically 
symmetric one-dimensional code. Using this code we can corroborate with high numerical accuracy 
that the constraint violations in the formulation of the equations that we use do not leak to the 
outside, as expected from our analytical analysis. 

The one-dimensional numerical studies also reveal an interesting property of turducken evolu- 
tions: even though the stuffing procedure initially introduces large constraint violations inside the 
black holes, these violations quickly decay to very small values as the evolution proceeds. This 
occurs because the shift vector quickly moves the coordinate grid points away from the future do- 
main of dependence of the turduckened region, while the constraint-violating modes are confined 
to the inside of the black hole. The numerical data then relax to a portion of the stationary 1 + log 
"trumpet slice" of the black hole 11, j3, 13]. This is the same end state as obtained with puncture 



evolution. 

We also use one-dimensional simulations to investigate the possibility of superluminal gauge 
modes. We find that gauge modes are in fact superluminal and propagate from the interior to the 
exterior of the black hole. In particular, the smoothing procedure affects the coordinate conditions 
outside the black hole. However, we find that the differences in gauge that arise from different types 
of smoothing quickly decay in time. As already mentioned, we find that the turducken solution 
approaches a portion of the trumpet slice, regardless of the type of smoothing. 

Given that it has already been shown in [Hj^j that the turduckening procedure works in practice 
for binary black hole evolutions, we next analyze in detail several aspects of single black hole 
evolutions. In Sec. |V] we present results from three-dimensional evolutions of a single distorted 
rotating black hole. This data is obtained by applying the smoothing procedure to puncture 
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initial data with Bowen-York extrinsic curvature. We compare in detail the gauge conditions and 
extracted waveforms produced in calculations with turduckening regions of different sizes as well 
as a pure puncture evolution. We also show that we have fourth order convergence in the extracted 
waveforms. 

We find, in agreement with the one-dimensional results, that superluminal gauge modes are 
able to propagate to the outside of the horizon. However, if the turduckening region is sufficiently 
small, the effect of these gauge modes decreases with radius outside of the black hole, and becomes 
small enough that for practical purposes it can be disregarded in 3D simulations. 

Comparing the waveforms from turduckening and pure puncture runs, we find that the dif- 
ferences are very small and that most of them converge to zero with increasing resolution. The 
remaining differences are caused by the differences in gauge at the finite detector locations and are 
found to become smaller with increasing detector location. 

Section IVII contains some final remarks. 



II. FORMULATION OF THE EQUATIONS, CONSTRAINT PROPAGATION, 
HYPERBOLICITY, AND CHARACTERISTIC SPEEDS 

In this section we first give the explicit form of the evolution and constraint equations used 
in our Cauchy formulation of Einstein's field equations. It is a special case of the family of for- 



mulations analyzed in 14J. Next, we summarize the conditions under which this formulation is 



hyperbolic and give the characteristic speeds. Finally, we extend the analysis performed in 14] 
by deriving the constraint propagation system, describing the propagation of constraint violations, 
and analyzing its hyperbolic structure. In particular, we give necessary conditions for this system 
to be symmetric hyperbolic and possess no superluminal speeds. We then prove that under these 
conditions constraint violations inside a black hole which are present in the turducken approach 
cannot propagate to the domain of outer communication. 

A. Formulation of the equations 

As mentioned in the introduction, we consider a BSSN-type formulation of Einstein's equations 
where the three metric 7^ and the extrinsic curvature are decomposed according to 

Hj = e 4 %j , (1) 
Kij = e 4 * (% + %K] . (2) 



Here, the conformal factor e 2 ^ is chosen such the conformal metric 7^ has unit determinant, and 
K = 7 JJ Kij and Aij are the trace and the trace-less part, respectively, of the conformally rescaled 
extrinsic curvature. The 3 + 1 decomposition of Einstein's equations alon g w ith suitable gauge 



conditions for lapse (a) and shift (/?*) yields the following evolution system 14J 1 



1 There are two sign errors in Eqs. (5) and (6) of Ref. [3- The first is in front of the second term of Eq. (5) and the 
second in front of the fourth term in Eq. (6). Since these errors only affect lower order terms they do not affect 
the results in in any way. We thank Dae-Il Choi for pointing out these errors to us. 
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where we have introduced the operator do = dt — 3 dj- Here, all quantities with a tilde refer to the 
conformal three metric 7™ and the latter is used in order to raise and lower indices. In particular, 
Di and T\j refer to the covariant derivative and the Christoffel symbols, respectively, with respect 
to jij. The expression [• • - ] TF denotes the trace-less part (with respect to the metric 7^) of the 
expression inside the parentheses, and 
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The gauge conditions imposed on the lapse, Eq. ([3]), is a generalization of the Bona-Masso 



condition [15[ where / (a, 1 



is a smooth and strictly positive function and Ko(x fl ) is an arbitrary 



smooth function. The conditions imposed on the shift in Eqs. (|5l6p is a generalization of the 



hyperbolic Gamma driver [16|] condition where G(a 



and H(a, < 



are smooth, strictly 



positive functions, and rf(B 3 , a, x^) is a smooth vector-valued function. The term doV 1 in Eq. (J6j) 
is set equal to the right-hand side of Eq. (fTOj) . Note that we use the operator do (as opposed to 
dt) in these gauge conditions; not only does this simplify the analysis of the principal part of the 
evolution equations, it also results in stable binary black hole evolutions for moving punctures 17] 
and the turducken approach [l| . 



Finally, the parameter m which was introduced in 18[] , controls how the momentum constraint is 
added to the evolution equations for the variable T l . The standard choice in numerical simulations 
is m = 1 which eliminates the divergence of A 3 in Eq. (|10p . However, we find it instructive not 
to fix m = 1 in this article. The source terms S, Si 3 and S l are defined in terms of the four Ricci 
tensor, R Y/ , and the constraint variables 



H ee \ + K 2 



K l3 K, 



Mi 



D 3 A i3 - ^D t K + QAyb 3 ^ 

r + djf 3 , 



(13) 
(14) 
(15) 
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as 



S = ~f ij R, 



(4) 



2H, 



S. 



i.i 
S' 



TF 



2amf j M j -d C l T . 



(16) 
(17) 
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In vacuum, the evolution equations consist of Eqs. pl fTUI) with 5 = 0, Sij = 0, S l = 0. In order to 



obtain a solution to Einstein's vacuum field equations, one also has to solve the constraints H = 0, 
Mi = and Cp = 0. Below, we show that for m = 1 it is sufficient to solve these constraints on 
an initial Cauchy surface in the region exterior to black holes. The constraint propagation system 
then guarantees that these constraints hold at all events which are future to the initial surface and 
outside the black hole regions, provided suitable boundary conditions are specified at the outer 
boundary of the computational domain. 



B. Hyperbolicity and characteristic speeds for the main system 



The evolution system (pl llOp is first order in time and mixed first/second order in space. There 
exist at least three different methods for analyzing hyperbolicity (that is, the well-posedness of 
the Cauchy formulation) for such systems. The first method consists in reducing the system to 
fully first order by introducing extra variables (and constraints) and to show that the resulting 
first order system is strongly or symmetric hyperbolic (see [3] for definitions). The hyperbolicity 
of the BSSN equations with a fixed shift and a densitized lapse or a Bona-Masso type condition 
using this method has been established in Refs. [2(3] and 14]. The second method which was 
developed in Refs. 
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221 ] is also based on a first order system. However, the reduction makes use 
of pseudo-differential operators. This has the advantage of not introducing any new constraints. 
Unlike the first method, this results in a unique first order system. The hyperbolicity of the BSSN 
equations with a Bona-Masso type condition and a hyperbolic Gamma driver type condition was 
shown in 14] using this method. Finally, the third method which has been proposed in Ref. [23| 
and applied to BSSN in Ref. 24J consists in finding an energy norm which, in the limit of frozen 
coefficients, is conserved. This method has been shown [25| to be equivalent to obtaining a first 
order symmetric hyperbolic reduction with the first method. 

Based on the second method, the following characteristic speeds with respect to normal observers 
for the evolution system ([3l ll0p were found in (l4T [: 0, ±/ii, ±/i2,±/U3, ±/M, ±^5 an d ±/i6 where 



Ml 



77, 



/'2 



4m — 1 



A*3 = Vm 



/i 4 = 1, fJ-5 = VGH 



/'6 



AGH 



(19) 
(20) 



When considering high-frequency perturbations of smooth solutions it is possible to classify the 
characteristic fields as gauge fields, constraint- violating fields and gravitational radiation [261. 1271]. 
According to this classification, the fields propagating with speeds fix, /is and correspond to 
gauge modes, the fields propagating with speeds [12 and 113 to constraint-violating modes and the 
fields propagating with speeds /i4 to gravitational radiation. As shown below, this statement can 
be strengthened by noticing that 0, /i2 and ^3 are the characteristic speeds of the constraint propa- 
gation system. In fact, it can be shown [28[] under quite general assumptions that the characteristic 
speeds of the constraint propagation system are a subset of the speeds of the main evolution system. 

In 14j] the following necessary conditions for strong hyperbolicity are given: f > 0, m > 1/4 and 
GH > or / > 0, m > 1/4 and G = H = 0. (Notice that for G = H = the evolution equation 



for the shift, Eq. ([5]), decouples from the remaining system.) If, in addition, the parameter m and 
the functions /, G and H can be chosen such that the functions 

AGH 6(m - 1) 2(m - l)GH 



3/ - AGH ' 4m - 1 - 4G# m-GH 
have smooth limits at 3/ = AGH, Am = 1 + GH and m = GH, then strong hyperbolicity is 



guaranteed [14J. For the standard choice m = 1 it is sufficient to verify that / > 0, Cff > and 
that the function AGH/(3f - AGH) has a smooth limit at 3/ = AGH. 

In the three-dimensional simulations below, we fix the functions /, Kq and G, H and rf as 
follows. We choose the 1 + log condition 

/ = - , ^o = 0, 



a 



and the Gamma-driver shift condition 

3 

4a 2 



with j] = 1/2. In this case, fi± = \j2ja and strong hyperbolicity is guaranteed if the function 
2ae -4 ^ — 1 does not cross zero. In our initial slices a — > 1 and <fi — > in the asymptotic region 
while near black holes a > is small (a ~ 0.3 at the horizon) and (j) positive. Therefore there 
already exists a two-surface where the condition 2ae~^ — 1 7^ is violated in the initial data. On 
the other hand, since this surface is a set of zero measure in the computational domain there is 
hope that the violation of our sufficient conditions at this surface might still result in a well posed 
Cauchy problem. The numerical simulations below show no apparent sign of instability. 



C. Hyperbolicity and characteristic speeds of the constraint propagation system 

Next, we derive the constraint propagation system which describes the propagation of constraint 
violations. We prove that for 1/4 < m < 1 constraint violations inside a black hole region cannot 
propagate to the outside. 

A convenient way of finding the constraint propagation system is to perform a 3 + 1 decompo- 
sition of the contracted Bianchi identities, 2V^R^} — \/ u R^ = 0, where one sets the quantities 
S, Sij and S l defined in Eqs. (|16H18p to zero. Taking into account the definitions of the constraint 
variables H , Mj and Cp defined in Eqs. (|13M15p one finds that they obey the linear evolution system 

d oH = _I D i(o? M )) - ae-^AV^djC^ + ^ KH, (21) 
a 3 



a 3 



TF 



doMj = —D 3 {a- 2 H) + aKMj + + D % \a [% { id j} C^ J , (22) 

d Cl = 2amfiMj . (23) 

In order to analyze this system, which is mixed first/second order in space, we use the first method 
described in Sec. Ill Bl and reduce it to a first order symmetric hyperbolic system. This allows us to 
establish the causal propagation of the constraints via a standard energy inequality. Introducing 
the additional constraint variable Zi k = 3jCp, Z\j = Z^^i-, Eqs. (|21H23p can be rewritten as the 
following first order linear system: 
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— D 3 (a Ms) 
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2a 



ae-^A ij Zij + — KH, 



■Dj(a- 2 H) + aKMj + M i d j p i + D i (aZ {ij) ) TF - aae~^ d k Zj k - djZ k 



2amf j Mj . 

2mdi (aMj) - 2amj kl (da jk )Mi 
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+ Z lk dj(3 k + Z kj di[3 k + Z k %3 k (5 l - ^Zijd k f3 k . 



(24) 

(25) 
(26) 

(27) 



Here, we have included in the right-hand side of Eq. (|25j) the term d k Zj k — djZ k k with an arbitrary 
factor a. Since Z k = diC k , this term is identically zero. However, as we will see now, its addition 
allows greater flexibility in obtaining a symmetric hyperbolic system. The system (1241127P has the 
form 



d C = a [A(u)%C + B(u)C] , (28) 

where C are the constraint variables, u = (a, f3 l , <f>, K,^j, Aij) are the main variables, and A 1 , 
i = 1,2,3, and B are matrix-valued functions of u. Decomposing Zy = Zuj) + Z\ij] +"fijZ/3 
into its trace-free symmetric part, Znj\, its antisymmetric part, Zuji, and its trace, Z = 7 4J Zij = 
e~^Z k k , and representing C in terms of the variables C = (C^,S\ := 2mH + Z, S2 '■= H + 
2aZ,Mj,Z(ij-),Zuj]), the principal symbol A(n) = A(u) l rii is given by 



A(n) 



/ ci \ 1 \ 

S x 
S 2 (4,ma - l)n j Mj 

\njS2 + (1 - a)n % Z^j) + an l Z[ij] 
2m(n (i M i) ) TF 
V%])7 V 2mn [i M j] 



Mj 



(29) 



Here n l = ^rij and re» is normalized such that run 1 = 1. This system is symmetric hyperbolic if 
and only if the following inequalities hold: 

4ma - 1 > 0, 2m(l - a) > 0, 2ma > 0, 

which is equivalent to the two conditions m > 1/4 and l/(4m) < a < 1. Therefore, as long 
as m > 1/4 (which is also a necessary condition for the main evolution system to be strongly 
hyperbolic, see Sec. Ill B| ) we can choose a between l/(4m) and 1 and obtain a symmetric hyperbolic 
constraint propagation system. For the standard choice m = 1, for instance, we can choose a = 1/2 
which is the case considered in [l||. A symmetrizer H = H T is given by 

C T HC = j l3 C l r Cl + Si + 



Si + Y J MiMi 

3(4mo- - 1) 2 J 



a 



^ l z m z ikl) + 



2m ' ' WJ ^ L > 2m 
H is positive definite and satisfies HA(n) = A(n) T H. The symmetrizer allows us to obtain an 
energy-type estimate 2 for the constraint variables C. For this, define the four-current 

M 2a 2 



2 Such estimates are a standard technique in the theory of hyperbolic partial differential equations. In particular, 
they allow one to prove uniqueness and continuous dependence on the data and to establish the principle of finite 
propagation speed. For references, see for instance [28l.l29||. 
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By virtue of Eq. (|28p the current satisfies the conservation law 

= dt J* + dif = C T LC, (30) 
with 2L = HB + (HB) T + d t (a^H) - d { (HA ! + a" 1 /?^) . Next, let ft T = [j E t be a tubular 

0<i<T 

region obtained by piling up open subsets T, t of t = const hypersurfaces. This region is enclosed 
by the space- like hypersurfaces So, £t and the surface T := |J <9£t, which is assumed to be 

0<t<T 

smooth. Integrating (|30p over $7^ and using Gauss' theorem in M. with the Euclidean metric, one 
obtains 

J J t d 3 x = j J t d 3 x + J C T LCd A x - J J%d5, 

where is the outward unit one-form normal to T and dS the volume element on that surface. If 
the boundary term J^e^ is positive or zero, one obtains the estimate 

J J t d 3 x < J J t d 3 x + J C T LCd A x 



C 1 LCd 3 x dt 




J f d 3 x I dt, 



where b is a constant and where we have used the positivity of J* = C T HC/(2a) in the last step. 
By Gronwalhs lemma, one obtains the inequality 



J t d 6 x < e bt J rd 6 x, 0<t<T. (31) 

Since 2aJ* = C r HC is positive definite, this then implies that C = everywhere on Q,t if C = 
on Eo which shows that it is sufficient to solve the constraints C = on the initial slice So- In 
view of numerical applications, however, the constraints are not exactly satisfied on So- Instead, 
numerical errors introduced by solving the constraint equations on a finite grid may be modeled 
by a sequence C n of initial constraint violations which converges to zero as the resolution goes to 
infinity. In this case, the estimate (I3ip shows that for each fixed t € [0, T] the L 2 -norm of the 
constraint variables C n converge to zero on 

In order to analyze the conditions under which the boundary term is nonnegative it is convenient 
to expand the outward normal vector as e^dx^ = N[a adt + rii(dx l + /3 l dt)] where rtj is normalized 
such that ^riifij = 1 and N > is a normalization factor. We set N = 1 in the following since we 
are only interested in the sign of the boundary term. Notice that \a\ < 1 if T is time-like, \a\ > 1 
if T is space-like and \a\ = 1 if T is a null surface. With this notation, the boundary term is equal 
to 

J^e M = \p T [aH - HAX] C. 

Therefore, the condition for this boundary term to be positive or zero is that all the eigenvalues of 
the matrix A*n, are smaller than or equal to a. If T is a future event horizon, then a = 1 and this 
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condition means that all of the eigenvalues must be smaller than or equal to one. For the symbol 
given in Eq. ([29]) these eigenvalues are the characteristic speeds (with respect to normal observers) 
and are 



0, ±fj,2 = ±d — - — , ±/i 3 = ± v / m. 

In particular, there are no super luminal speeds if 1/4 < m < 1, and in this case no constraint 
violations can propagate out of a black hole. For the standard choice m = 1 this condition is 
satisfied. However, we also see that choosing m > 1 in the evolution equation for T l , Eq. (I10D . 
yields superluminal constraint speeds in which case constraint violations inside a black hole can 
affect the exterior region. 

Finally, we would like to point out that if the computational domain contains time-like bound- 
aries, then \a\ < 1 and the sign of the boundary term J^e^ is not automatically positive or zero. 
In this case, boundary conditions need to be specified such that this term is nonnegative and such 
that a well posed Cauchy problem is obtained. 



D. Spherical symmetry 



The BSSN equations can be specialized to spherical symmetry as described in Ref. 30j. The 
first step is to remove the restriction 7 = 1 on the determinant of the conformal metric and replace 



it with an evolution equation for 7 31)]. In this paper we use the "Eulerian evolution" defined by 
dt 7 = 27 Difi 1 . The conformal connection functions are defined in terms of the conformal Christoffel 
symbols by T l = J^ k Tj k - 

The reduction to spherical symmetry is achieved by writing the conformal metric as ds 2 = 
^rrdr 2 + ; ygQd 2 Q, where d 2 Sl is the metric for the unit two-sphere. The independent component 
of the trace-free part of the extrinsic curvature is A rr , and the independent component of the 
conformal connection functions is r r . 

In one dimension as in three, we use 1 + log slicing and the Gamma-driver shift condition 
(although here we use r/ = 1 for the damping parameter). For the main evolution system in 
spherical symmetry the characteristic speeds are 0, db/zi, =t/i2 and ±^5. Strong hyperbolicity is 
guaranteed for / > 0, m > 1/4 and GH > if the following conditions hold: / / GH and 
3GH + 1 + 4m. 

The constraint evolution system can be obtained by spherical reduction of the system (1211123^ . 
or by direct calculation from the ID equations of motion Let = <9 r Cf and define the 

vector of constraints by C = (H, M r ,C£, Z^) T . The constraint evolution equations have the form 
8qC = a[A r 3 r C + BC] where A r and B are functions of the BSSN variables. The principal symbol 
is given by 



/ -2e 

1/6 



\ 2m/7 rr 



7r 





2e~ 4< ? 







/3 



\ 



(32) 



/ 



The characteristic fields are Cf, mH + e 



-^Z;, and H ± 2^3(4™ - l)/7 rr e _2 *M r + Ae'^Z^ with 



proper speeds 0, 0, and ±//2 = ±y (4m — l)/3, respectively. This system is symmetric hyperbolic 
as long as 4m > 1. For the case of primary interest, m = 1, the characteristic fields propagate 
along the normal and along the light cone. 

A symmetrizer can be constructed from the squares of the characteristic fields: 
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C T UC = (mH + e~^Z r r f + (Cf ) 2 + (# + 2^3(4™ - l)/ 7rr e- 2 *M r + Ae~ 4 * Z r r f 
+{H - 2V3(4m - l)/% r e-^M r + 4e~^Z r r ) 2 
= (m 2 + 2)H 2 + (Cf ) 2 + 24(4m - l) e -^(M r ) 2 /7 rr 

+33e^(Z r r ) 2 + 2e~ 4 *(m + 8)FZ r r . (33) 

The spacetime current, denned by 

J»d u = — C T UC d - -C T UA r C d r , (34) 

" 2a 2 ' v ; 

satisfies the conservation law 

= d t J l + d r J r = C T LC (35) 

where 2L = HB + B T H + d t (H/a) - <9 r (H(A r + (i r /a)). As in the 3D case, we can now show that 
if J M e^ is non-negative at the boundaries, then 

I J t dr<e bt J J l dr , < t < T (36) 
J Tit ■'So 

for some constant b. It follows that the constraints will vanish on £j if they vanish on the initial 
hypersurface So- 

The assumption that J^e^ is non-negative at the boundaries for m = 1 can be seen to hold at the 
black hole horizon by following the same reasoning as in the three-dimensional case. We expand the 
Euclidean normal one form as a linear combination of the unit normal to the t = const surfaces 
and the unit normal rii to the two-dimensional boundary within the spacelike hyper surf aces. As in 
the three-dimensional case, we have e^dx^ = N[aadt + n r (dx r + f3 r dt)] where \a\ = 1 characterizes 
a null surface. By dropping the positive constant N, we find J^e^ = C T (aH — HA r n r )C/2. This 
shows that for a black hole horizon J^e^ is positive if the eigenvalues of A r n r are less than or equal 
to one. This is indeed the case for m = 1, since the constraint propagation system has eigenvalues 
and ±^(4™ - l)/3. 



III. CODE DESCRIPTIONS 



We use two different codes for the simulations presented in this paper. One of them is 
McLachlan, a three-dimensional adaptive mesh refinement code, which uses the BSSN system 
of equations, as described in If], 32], with the gauge conditions described in section [TT1 above. See 
section [TT] for the exact form of these equations. The McLachlan code is a Cactus thorn which is 
entirely generated by Kranc H, 34, HH] directly from equations and differencing stencils specified 



in Mathematica format. McLachlan uses the Cactus framework [36l . l37| | and the Carpet mesh re- 
finement driver 



39(. The evolution code is fully fourth order accurate in time and space. We 



use fifth order spatial interpolation at mesh refinement boundaries, using buffer zones as described 



m 



to ensure stability and convergence at mesh refinement boundaries, and using tapered grids 
as described in [40j] to avoid the Berger-Oliger time interpolation. Our finite differencing opera- 
tors are the standard centered fourth order accurate first and second finite differencing operators, 
except for the advection terms which are upwinded (also fourth order). We use a fourth order 
Runge-Kutta time integrator and add fifth order Kreiss-Oliger dissipation [4l( to the right hand 
sides. We apply standard radiation (Sommerfeld) boundary conditions (as described in [3]) to 
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all components of the evolved fields. These boundary conditions are neither 4th order convergent, 
nor constraint preserving so non-convergent constraint violations will propagate inwards from the 
outer boundary. We therefore expect our code, in the limit of infinite resolution, to be fully fourth 
order accurate only in the region that is causally disconnected from the outer boundaries. We 
place our outer boundaries far enough out that they do not affect our wave extraction procedure. 

We center a stack of refined regions around the origin. Each region is a cube with a resolution 
half of the next coarser region. 

We also use a one-dimensional code to analyze several issues related to the turducken technique 



in the context of a spherically symmetric black hole. This is the code used in [12|] and described 
in detail as the "Eulerian case" in [3(j. The ID code uses a uniform radial grid with nodes at 
coordinate radii rj = (j — 1/2) Ar, where j = 1, 2, etc. Fourth order finite differencing is used 
for spatial derivatives, and fourth order Runge-Kutta is used for the time update. No boundary 
conditions are imposed at the origin for any of the variables except the shift vector component f3 r . 
That is, for all variables except (3 r , the finite difference stencil is shifted toward positive r near 
the origin so that no guard cells are needed. For f3 r we impose the boundary condition d r (3 r = 
at r = by using a fourth order, one sided finite difference representation of d r (3 r . The resulting 
guard cell value is 

/r(o) = 1 [i7/r (i) + 9/r(2) - 5/r(3) + f{a)\ (37) 

where the numbers in parentheses label grid points. In the evolution code, spatial derivatives of (5 r 
are computed by shifting the finite difference stencil toward positive r near the origin. In this way 
only the single guard cell value (3 r (0) is needed. Numerical experiments show that the condition 
(I37p (or a similar one) is needed for stability. 

The one-dimensional code uses the variable x — rather than <fi. The CFL (Courant- 

Friedrichs-Lewy) factor is 0.25, and unless otherwise stated the resolution for all the one- 
dimensional runs discussed here is Ar = M/200. 



IV. SINGLE BLACK HOLE EVOLUTIONS AND THE END STATE 

In this section we investigate the effects of black hole smoothing on the evolution of a spherically 
symmetric single black hole. 

We start with the (Schwarzschild) isotropic black hole data 

7 rr = 1 (38a) 

lee = r 2 (38b) 

= l + M/(2r) (38c) 

f r = -2/r (38d) 

along with K = and A rr = 0. For all one-dimensional simulations, we use the 1 + log slicing and 
Gamma driver shift conditions as described in Sec. IIIB1 The initial values for the gauge variables 
are a = 1 and (3 r = B r = 0. Unless otherwise stated, all simulations use m = 1. 

The black hole interior is turduckened by making the replacement r — > f(r) in the data (|38p for 
< r < rt, where 



f(r) = ro - (10r - 4r t )r 2 /V t 2 + (20r - 6r t )r 3 /rf - (15r - 4r t )r 4 /rf + (4r - r t )r 5 /r% . (39) 
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The function r has the properties f(0) = ro, f'(0) = 0, f(r t ) = r t , f'(r t ) = 1, f"(r t ) = 0, and 
f'"{rt) = 0. Thus, with this turduckening, we extend r inside r^ in such a way that the initial 
data are C 3 at rt and nonsingular at r = 0. The lack of smoothness at rt generates an error 
in any centered fourth order finite difference derivative whose stencil extends across rt- For first 
derivatives the error is 0(Ar 3 ), and for second derivatives the error is 0(Ar 2 ). Note that initially 
the horizon is located at coordinate radius r = M/2. As long as rt < M/2, the turduckening lies 
entirely inside the black hole. 

The turducken runs considered below use either ro = 0.1M, rt = 0.4M or ro = 0.05M, n = 
0.45M. These two types of smoothing are referred to as case TA and case TB, respectively. We 
compare these results to results obtained with puncture data, denoted P, which is equivalent to no 
turduckening (puncture evolution). 

A. Behavior of constraint violations 

Initially the constraint violations are restricted to the black hole interior, and our analysis in 
Section [II] shows that for m < 1 they should stay there. Thus, we expect the constraints to hold 
everywhere and at all times in the black hole exterior. 

One obvious way to test this is to monitor the numerical constraints as functions of time, and 
confirm that the violations introduced by the turduckening do not leak out of the black hole. This 
is indeed the case. For smoothing types TA and TB, the initial Hamiltonian constraint violation 
inside the black hole is ~ 10 _1 /M 2 , while (at resolution Ar = A//200) the initial constraint 
violation outside the black hole is ~ 1CP 9 /M 2 . The constraint violation outside the black hole 
remains ~ 10~ 9 /M 2 throughout the evolution. 

Moreover, we find that the region of constraint violation quickly shrinks relative to the nu- 
merical grid, and the constraints quickly loose memory of the turduckening. This comes about 
because the grid points surrounding the origin acquire a radially outward velocity that becomes 
superluminal within a time of a few M. The curves labeled "coord" in Figs. Q] show the proper 
speed of the coordinates in the radial direction with respect to observers at rest in the spacelike 
hypersurfaces. The speed is plotted as a function of proper distance from the black hole horizon, 
with the convention that positive values are outside the black hole and negative values are inside 
the black hole. What we see from these figures is that the coordinate grid inside the black hole 
moves faster than the speed of light (c = 1) in the radially outward direction. The region of con- 
straint violation moves causally, with speeds and ±1. Thus, the coordinates soon move outside 
the future light cone of the stuffed region, and into the forward light cone of the initial data that 
satisfies the constraints. The graphs in Fig. Q] were taken from simulations with smoothing type 
TA. The graphs obtained with smoothing TB are nearly identical. 

In Figs. [2] we plot the Hamiltonian constraint as a function of coordinate radius for the initial 
data t = and for the later times t = AM, 8M, and 12M. The data for turduckening types TA 
and TB are shown, as well as for puncture data P. The constraint violation does not propagate 
beyond the black hole horizon. The value of the Hamiltonian constraint violation inside the black 
hole drops from ~ 0.1 to ~ 10~ 9 by t = 12M. Beyond about t = 12M, the constraint data for TA 
and TB are everywhere (that is, also inside the black hole) indistinguishable from each other, and 
indistinguishable from the results obtained with puncture data P. 

The success of the turduckening procedure depends on the constraint violating modes propa- 
gating with speeds less than or equal to one. Our analysis in Section [II] predicts that this condition 
is met for the BSSN family of evolution equations ([3j)- (ll0p if and only if 1/4 < m < 1. In Fig. [3] 
we plot the numerical Hamiltonian constraint as a function of coordinate radius r for m = 1.25. 
One can clearly see that, as predicted by the theory, in this case the constraint violation does 
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Figure 1: Proper speeds as a function of proper distance from the horizon for the turduckening case TA. 
The curve labeled 'coord' is the proper speed of the coordinate system with respect to the normal observers, 
^^f rr e 2 ^P r '/a. The curve labeled "mode 1" is the proper speed [i\ = \/2/a. The curve labeled "mode 5" is 
the proper speed ^5 = v / 3e 2( ^/(2a). The horizontal line is light speed. 



propagate from the interior to the exterior of the black hole. This illustrates the fact that for a 
given formulation of the Einstein equations one cannot simply assume that the constraints will 
propagate with non-superluminal speeds. 



B. Behavior of the coordinates 

The relatively large value of the radial component of the shift vector moves the grid points 
beyond the region of constraint violation within a time of a few M. However, the grid points 
do not move beyond the influence of the turduckening. Recall that the main evolution system in 
spherical symmetry is strongly hyperbolic with characteristic speeds 0, ±/ii = iy^/a, ±/i2 = ±1, 
and ±/X5 = ±\^3e 2c l > / (2a) . The modes with speeds \i\ and ^5 can become superluminal. Figure Q] 
shows the speeds /ii and fi§ at times t = 3M and t = 6M, along with the speed of the coordinate 
grid relative to the normal observers. Both speeds [i\ and are larger than the coordinate speed 
at the black hole horizon (the origin of proper distance). As discussed in Sec. IIIB| in some sense 
the modes corresponding to fi\ and ^5 can be associated with gauge freedom. We therefore expect 
that the turduckening process can affect gauge conditions outside the black hole. Although the 
full Einstein equations are satisfied in the black hole exterior, independent of the smoothing, the 
slicing and coordinate conditions outside the black hole can depend on the details of the stuffing. 

Figure H] shows the differences between lapse functions versus areal radius R, at four different 
times. For the curve TA-P, we compute the difference between the lapse function for smoothing 
TA and the lapse function obtained from puncture evolution (no smoothing). These calculations 
require some explanation. The areal radius is not a monotonic function of the coordinate radius; 
rather, the areal radius has a minimum at the black hole throat. Therefore for each of the runs TA 
and P we use only the data for which R is an increasing function of coordinate radius. We then 
consider the overlap region in which R is increasing for both data sets TA and P. The difference 
in lapse values is computed by interpolating this data onto a common grid that covers the overlap 
region. The differences TB— P are computed in the same way, but with data from smoothing type 
TB. Note that the overlap regions extend inside the black hole horizon, which has areal radius 
R = 2M. 
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distinguished from the horizontal axis. shows some small differences that disappear shortly 

after t = 12M. 



Figure 2: Hamiltonian constraint versus coordinate radius r at various times, in units of M. TA and TB 
denote two types of turduckening, and P is puncture data. The vertical line shows the location of the black 
hole horizon. 



The effect that the turduckening has on the slicing is relatively small when compared to the 
nominal value of 1 for the lapse, but is clearly seen in Fig. |H This effect begins inside the black 
hole and propagates super luminally to the outside, where it continues to spread radially outward. 

The turduckening's effect on the slicing condition fades with time. As the evolution proceeds, 
the data relax to a portion of a stationary 1 + log slice of a Schwarzschild black hole, independent 
of the initial stuffing details. This stationary 1 + log slice has a "trumpet" geometry [l3| ■ It is the 



same final slice obtained with puncture evolution [12, |42j, |43| . In Fig. [5] we graph the areal radius 
as a function of proper distance from the horizon (with the convention that positive distances are 
outside the black hole, negative distances are inside). These plots show the data for smoothing 
types TA and TB, and for no smoothing P, at early (t = 0.25M) and late (t = 50M) times. The 
stationary 1 + log slice of Schwarzschild is shown as the curve labeled S. Initially the R versus 
proper distance relation shows a strong dependence on smoothing. By t = 50M all of the data 
TA, TB, and P have relaxed to approximate a portion of the trumpet slice S. At this time, and 
with resolution Ar = M /200, the numerical slices all end at proper distance ~ — 6.45M. Close 
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(a)Time t = AM. With m = 1.25, the constraint 
violation is beginning to spread to the black hole 
exterior. 



(b)Time t = 8M. With m = 1.25 a pulse of constraint 
violating data has moved beyond the black hole and 
continues to propagate outward. 



Figure 3: Hamiltonian constraint versus coordinate radius r for m = 1 and m = 1.25, where m is the 
parameter that controls the mixing of the momentum constraint in the equation of motion for r r . The 
smoothing type TB is used in both cases. The vertical line shows the location of the black hole horizon. As 
predicted by the theory, the constraint violations inside the black hole do propagate to the outside if m > 1. 



inspection of the data shows that near the end of the numerical slice, the areal radius for the cases 
TA, TB, and P agree with one another to more than seven decimal places; at proper distance 
— 6M, we find a ~ 1.32013. The areal radius for a stationary 1 + log slice at proper distance — 6M 
is a w 1.32018. 

Although the stuffing can affect the slicing beyond the black hole horizon, it does not always do 
so. For the type of simulations here considered, it appears that any stuffing that is initially inside 
a coordinate radius of about 0.2M remains causally disconnected from the black hole exterior. In 
Fig. Owe show the difference between lapse functions for stuffing types TD, TE and puncture data 
P. The cases TD and TE use turduckening radii rt = 0.2M and r% = 0.3M, respectively. Both 
cases use ro = 0.02M. The figure shows the common logarithm of the lapse differences at t = 3M 
and t = 6M, for low (Ar = M/100) and high (Ar = M/200) resolutions. 

The top two subfigures show that the lapse difference TD— P converges to zero with increasing 
resolution. On the log plots the difference between low and high resolution curves is log(16) ~ 1.2. 
These results show that with its small turduckening radius, the data for stuffing type TD are 
indistinguishable from puncture data in the black hole exterior. The bottom two subfigures show 
the lapse difference TE— P. In this case the difference converges to a nonzero value which spreads 
from the black hole interior to its exterior. 

Let us refer once again to Fig. [H which shows the gauge propagation speeds fj,i, along with 
the coordinate system speed. Observe that the speed \i\ = ^f] depends on the slicing condition 
([3]), while the speed /is = V GH depends on the coordinate shift conditions ([5]) and ([6]). Thus it is 
the speed /ii that we examine closely here. Figure Q] shows that the coordinate system moves faster 
than fi\ within a proper distance of ~ M of the puncture. This appears to be a common result, 
independent of the stuffing details. If the stuffing initially extends beyond this /ii-sphere, where 
the speed curve for m crosses the coordinate speed curve, then the stuffing's affect on the slicing 
can propagate radially outward through the black hole horizon. This is the case for data type TE. 
If the stuffing is initially contained entirely within the /xi-sphere, then the stuffing's affect on the 
slicing is lost as the coordinate grid quickly moves beyond the influence of mode fi±. This is the 
case for data type TD. 
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Figure 4: The difference between lapse functions plotted versus areal radius in units of M. The curve TA— P 
is the difference for turducken data with stuffing TA and puncture data. The curve TB— P is the difference 
for turducken data with stuffing TB and puncture data. 
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Figure 5: Areal radius R versus proper distance from the horizon. The horizontal line at R — 2 is the 
horizon. By t = 50M the three curves (TA, TB, and P) are indistinguishable (at this scale) from the 
stationary 1 + log slice S. 
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Areal radius Areal radius 

(c)Time t = 3M, turducken radius rt = 0.3. (d)Time t — 6M, turducken radius rt = 0.3. 

Figure 6: Common logarithm of the difference between lapse functions. The top two subfigures show the 
difference between case TD and puncture data for low (Ar = M/100) and high (Ar = M/200) resolutions. 
The bottom two subfigures show the difference between case TE and puncture data for low and high 
resolutions. (Note that the curves at time t = 3M appear to terminate at an areal radius of about 3.3 
(case TD) or 3.5 (case TE). This occurs because, beyond these values, the lapse difference is exactly 0.0 to 
thirteen or more decimal places. The logarithm is undefined for larger values of areal radius.) 

Note that the precise value of the turduckening radius where the influence of the gauge speed 
/ii is lost might depend on a number of details of the simulation. For example, black hole spin, 
the initial values of the lapse and shift, as well as the details of the stuffing profile might affect 
whether or not the gauge modes will be lost in the black hole interior. 

Finally, let us observe that the gauge speed /is is larger than the coordinate speed throughout 
the black hole interior. 3 Thus, we expect that the stuffing details will affect the shift in the black 
hole exterior, regardless of how small the stuffing region might be. This expectation is supported 
by the results shown in Fig. This figure shows the difference in coordinate radius for turducken 
data TD and puncture data P as a function of areal radius. As discussed above, in the limit of 
infinite resolution the slicing outside the black hole appears to be identical for turducken data TD 
and puncture data P. Thus the difference shown in Fig. [7J is due to differences in the coordinate 
grid. In particular, the radial shift in the coordinate grid in the black hole exterior depends on the 



3 This observation partially explains why the one-dimensional code [30| seems to require boundary conditions at the 
origin for the shift vector. 
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type of stuffing used. At time t = 2M the difference has just begun to cross the black hole horizon. 
The difference continues to move into the black hole exterior as the evolution proceeds. 
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Figure 7: The difference between coordinate radius for turducken data TD and puncture data, as a function 
of areal radius. Both low (Ar = M/100) and high (Ar = JI//200) resolution results are shown, although 
the curves are difficult to distinguish. The vertical line is the location of the black hole horizon. 



V. THREE-DIMENSIONAL BLACK HOLE EVOLUTIONS 
A. Turduckening procedure in 3D 

Numerical initial data with black holes may contain singular "puncture" points or excised 
regions. If the data are incomplete due to excision, the excised regions must be filled before 
the turducken method can be used to evolve the data. Even if the data are complete, as in the 
puncture case, it may be advantageous to replace the data near singular points with data that is 
more smooth. 

We experimented with various methods for turduckening the initial data in the black hole 
interior. We found that the details do not matter much in practice, as long as the fields are 
sufficiently smooth across the turduckening boundary and the spacetime remains unchanged within 
a layer inside the horizon whose width is at least a few times the width of the finite differencing 
stencil. Empirically, a width of 10 grid points suffices for us; we expect this number to depend on 
the particular differencing operators which are used. 

One rather simple method for turduckening is blending, which fills the excised region with 
arbitrary data, and then modifies some of the non-excised grid points to create a smooth match. 
This has the disadvantage that it may require quite a few non-excised grid points inside the 
horizon — one needs to have sufficiently many grid points to ensure the smooth match, plus the 
grid points which need to be left unmodified. 

Instead, we choose a method which leaves all given data unchanged and fills in all excised points 
in a smooth manner. In particular, we solve an elliptic equation of the form 

/ gn gn gn \ 

to fill the excised points of a quantity A, using standard centered derivatives everywhere and using 
the given non-excised data as boundary conditions. Here n is an even integer which controls the 
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Turduckening 


Resolution 




radius 


low med high 




n 


6A 4A 3A 


small 


0.10M 


11 17 23 


medium 


0.15M 


9 14 19 


large 


0.20M 


7 11 15 



Table I: Number of grid points between the turduckening region and radius of the initial apparent horizon 
(rhor = 0.3758M), for three turduckening radii r t and three grid resolutions dx. The resolution dx is given 
in multiples of A = 0.04M. If there are too few grid points between the horizon and the turduckening 
region, information about the turduckening procedure escapes out of the black hole, which must be avoided. 



smoothness of the resulting field A. If n = 2, A will be only C . Increasing re by 2 results in one 
additional derivative being continuous, so that for re = 4, A is C 1 , and for n = 6, the resulting A 
is C 2 . Since we need to take two derivatives of the metric, we choose n = 6 to keep all derivatives 
continuous. We will show later that using n = 2 still works, but leads to large errors. 

Equation (|40p is linear, and we employ a standard conjugate gradient method [4jy] to solve it. 
This numerical scheme is rather easy to implement, is robust, and it converges reasonably quickly 
at the resolutions within reach on current supercomputing hardware. The algorithm has been 
implemented in the Cactus thorn NoExcision that we have made freely available. 4 

Higher values of n lead to smoother initial data. It is also possible to choose a non-zero right 
hand side in (|40p . modifying the shape of the solution in the excised region. In this paper we do 
not take advantage of this additional freedom. Other turduckening procedures may also be possible 
and could be directly integrated into the initial data solvers. 



B. Distorted rotating black hole evolutions 

We now turn to evolutions of single distorted, rotating black holes and investigate the effects 
of black hole turduckening on the numerical solution and the extracted waveforms. We use single 



puncture initial data with a Bowen-York extrinsic curvature [45(. We use a puncture mass m p = 
0.751744 and a moderately large angular momentum parameter S = 0.7, resulting in a black hole 
with irreducible mass M irr « 0.925785, dimensionless spin a/M = S/M 2 = 0.7, and mass M = 1.0 
(see e.g. Eq. (27) in [46] for a definition of the horizon mass used here). The ADM mass of the 
spacetime is Madm ~ 1.00252. 

We chose puncture initial data to be able to compare turducken and puncture evolutions (see 
also Sec. HVp . We chose a non-zero angular momentum to arrive at a more interesting, non- 



spherically-symmetric case which also includes gravitational wave emission, as the conformally flat 



rotating single punctures do not represent the Kerr spacetime (4 71 ]. 

We perform simulations of four different initial data setups. The first is a pure puncture setup 
where the turduckening procedure is not applied. In the other three cases, the puncture data are 
modified inside of different turduckening radii: small (rt = 0.1M), medium (n = 0.15M), and large 
(rt = 0.2M), but kept unchanged everywhere else. The radius of the initial apparent horizon is 
r hor = 0.3758M. For all four cases we perform simulations at three different resolutions: low (dx = 
0.024M), medium (dx = 0.016M) and high (dx = 0.012M), as measured on the finest refinement 
level. TableUlists the approximate number of grid points between the horizon and the turduckening 
region in each of these cases. These runs were performed using the McLachlan code (see Sec. IIIip . 



4 The Cactus thorn can be obtained via the command 
svn checkout https : //svn. aei .mpg. de : /numrel/AEIThorns/NoExcision/trunk NoExcision 
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Figure 8: Convergence plot for the Hamiltonian constraint along the x axis at T = 115.2M for the small 
turduckening region rt = 0.1 M (left plot) and for the puncture run (right plot). The medium and high 
resolution curves have been scaled for 4th order convergence. The vertical line shows the location of the 
horizon. The plots cover the range over which the waveforms were extracted (up to R = 60A/) . At this time 
this region is still causally disconnected from the outer boundary. 



eight levels of mesh refinement, and the outer boundaries located at R = 256M. Refinement 



boundaries were placed at R = [128, 64, 16, 8, 4, 2, 1]M. Fifth order Kreiss-Oliger dissipation [41] 
was applied to all evolved variables. 
We choose the initial lapse profile as 

° init = l + 4/(2r)' (41) 

which corresponds to the average of the isotropic lapse and unit lapse. Except for the puncture 
case, the lapse is further modified in the turduckening region by smoothing it using Eq. (|40p with 
n = 6. The shift and its time derivative are initially set to zero. 



C. Numerical results 

1. Constraints 

Figure [8] shows the convergence behaviour for the Hamiltonian constraint for the small turduck- 
ening radius rt = 0.1M and the puncture run along the x axis at a late time T = 115. 2M (any 
late time could have been chosen, since the numerical spacetime has essentially become station- 
ary). The vertical line shows the location of the horizon at this time. Similar results hold for the 
momentum constraint. These results show that at late times the constraints converge to zero at a 
fourth order rate, except for a few grid points near the origin. Since the stationary solution we are 



approaching has a l/\/r singularity in the conformal factor 12j, |l3j, this is to be expected, as a 
neighborhood of the origin is under-resolved. Note also that the convergence plots are practically 
identical for the turducken and puncture simulations. These results are consistent with the ID 
results of Section IIV1 



2. Lapse 



The behavior of the lapse is also consistent with the spherically symmetric simulations of Sec- 
tion IIVI Figure [9] shows the lapse profile along the x axis at four different times (T = 0M, 
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Figure 9: Lapse profile along the x axis for the distorted single black hole evolutions at the times T = OM 
(top left), T = 3.84M (top right), T = 9.6M (bottom left), and T = 96M (bottom right) for the high 
resolution, comparing different turduckening radii. The vertical lines show the location of the apparent 
horizon. The insets in the bottom graphs enlarge a small region near the horizon. While the lapse profiles 
differ near the origin, they are very similar at the horizon, and their difference decreases with time. However, 
for the medium and large turduckening radii, a supcrluminal gauge mode is clearly escaping. On the other 
hand the differences between the lapse profiles of the puncture and small turduckening radius cases converge 
to zero with resolution. 



T = 3.84M, T = 9.6M, and T = 96M) for the high resolution case. The region with significant 
differences in the slicing is at all times safely contained within the horizon. However, as can be seen 
from the inset in the bottom left graph, there are real (but small) differences in the lapse function 
also outside the horizon at early times. These differences do not converge away with resolution. In 
other words, the details of the turduckening procedure result in real (but small) differences in the 
slicing outside of the horizon. With time, the differences become smaller, and are no longer visible 
in the bottom right graph at T = 96M. 



3. Waveforms: convergence with resolution 



We extract the t = 2, m = mode of the Weyl scalar ^4 on coordinate spheres at four different 
radii (R = 30M, R = 40 M, R = 50M, and R = 60M), choosing the commonly used hypersurface- 
adapted tetrad described e.g. in [48(]. Figure flQl shows convergence tests for the waveforms obtained 
with different turducken radii and with puncture evolutions, in all cases extracted at R = 30M. 

The top left graph shows the waveform with small turducken radius r$ = 0.1M at three different 
resolutions. The top right graph shows the differences between the different resolutions, scaled for 
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Figure 10: Effect of the turduckening radius on the i = 2,m = mode of gravitational waveforms 
extracted at R = 30 M . The top left graph compares all three resolutions, and the top right graph shows 
their differences, scaled for 4th order convergence, both for the small turduckening region r t = 0.1M. The 
bottom left graph shows the scaled differences for a pure puncture run, and the fact that it looks virtually 
identical indicates that r t = 0.1M is a good choice for these resolution. On the other hand, the bottom 
right graph shows the same for the large turduckening radius r = 0.2M where noticeable differences arc 
visible. We attribute this to the smaller number of grid points between the turduckening region and the 
initial apparent horizon (sec tabic [J), so this case is not yet in the convergent regime. 



fourth order convergence. The fact that the curves look so similar is a manifestation of clean fourth 
order convergence. 

The bottom left graph shows the scaled differences for a pure puncture run, also indicating 
clean fourth order convergence. The case rt = 0.15M (not shown here) behaves similarly. 

However, the scaled differences for rt = 0.2M (bottom right graph) looks noticeably different. 
We attribute this to the smaller number of grid points between the turduckening region and the 
initial apparent horizon (see tableU]), so that, in the low resolution case, there are too few gridpoints 
between the turduckening region and the horizon for the numerical waveforms to be effectively 
isolated from the stuffing. The curve for the difference between the medium and high resolution 
waveforms is very similar to the other cases, which indicates that it is only the low resolution run 
that is not in the convergent regime. We expect, though we have not yet confirmed this with 
numerical experiments, that a convergence test with a higher minimal resolution will also show a 
lack of convergence if the turduckening region is chosen so that rt is only ~ 7 (or fewer) gridpoints 
away from the horizon at the minimal resolution. Thus, as a rule of thumb, we expect that one 
should always choose the turduckening region so that r± is more than ~ 7 gridpoints away from 
the horizon. Finite differencing schemes that differ from the one used in our code would probably 
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Figure 11: Difference between the I = 2,m = mode of r^>4 for a puncture evolution and for three different 
turduckening radii, extracted at R = 30M (left-hand graphs) and R = 60M (right-hand graphs). See the 
main text for the discussion. 



require different limits. Experimentation would be necessary on a case-by-case basis. 



4- Waveforms: comparison among different solution methods 

Figure [TT1 shows, for the three resolutions, the differences between the puncture waveform and 
the turducken waveforms with turduckening radii rt = 0.1M, rj = 0.15M, and rt = 0.2M. In all 
cases two extraction radii are used: R = 30M and R = 60M. As can be seen from the figure, the 
difference between the puncture waveform and the rj = 0.1 one goes to zero as the resolution is 
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increased. This is not true for the rt = 0.15 and rt = 0.2 cases, where the difference, especially clear 
at R = 30M, does not converge to zero with increasing resolution. Differences in the slicing do 
appear to affect the extracted gravitational wave signal. However, at the larger extraction radius 
R = 60M, the non-convergent part of the waveform difference is significantly smaller. 

As already discussed, there are gauge modes that travel at superluminal speeds with our pa- 
rameter choices and differences in gauge are able to propagate beyond the horizon. This is true 
independently of whether the initial data are pure puncture or turduckened data. The fact that 
we only find differences between the puncture waveforms and turducken ones with larger rt does 
not mean that we consider the puncture waveforms and the turducken waveforms with small rt to 
be correct, and any deviation from them to be an error. As described above, we modify the initial 
lapse profile away from the puncture profile only inside of the turduckening region. We would 
expect that a pure puncture run with a different initial lapse profile (for example, constant lapse 
equal to one) would also yield significant differences in the lapse at the location of the detector. 
Instead, the conclusion is that since wave extraction is carried out at a finite radius, our resolutions 
are high enough so that the seemingly small differences in slicing lead to noticeable differences in 
waveforms. Furthermore, even though this is a rotating black hole, which cannot be compared 
directly to a Schwarzschild black hole, we believe that this is the same effect as seen in the ID 
spherically symmetric case (see Fig. [6] and its discussion in section UVBj) . Namely, differences in 
slicing get trapped inside the black hole for sufficiently small turducken radius (while the puncture 
evolution can be considered as the limit rt — > 0). Note also that the slicing differences (and conse- 
quently the waveform differences) are transient since both puncture and turducken runs approach 
the same trumpet slice at late times. 

What we have seen is that, through a proper convergence test, we can detect differences in the 
waveforms due to different slicings at fixed, finite extraction radii. However, it is far from clear that 
these differences are of any practical importance. For example, for each of our turduckening radii, 
the largest difference (at the highest resolution) between a turducken waveform and the puncture 
one is 6 x 10~ 6 . The amplitude of the wave is about 0.015, so the relative difference is about 
4 x 10~ 4 . 

Finally, we performed experiments at low (dx = 0.024M) and medium (dx = 0.016M) resolu- 
tions with a turduckening radius rt = 0.1, but with n = 2 in equation (|40p . so that the evolution 
fields are only C° at the boundary of the turduckening region. In Figure [12] the left graph shows 
the difference between the waveforms from runs with n = 2 and n = 6 while the right graph shows 
the difference between the n = 2 and puncture waveforms. These two graphs are almost identical 
and a comparison with the top left graph in Figure [TT1 explains why. As can be seen, the difference 
between the puncture waveforms and the rt = 0.1 waveforms is about 10 times larger for n = 2 
than n = 6. Note that in the rt = 0.1M case the turduckening region is far enough inside the 
black hole that we do not see gauge differences, so this must be caused by increased numerical 
noise coming from the less smooth data in the n = 2 case. 

VI. FINAL REMARKS 

In this paper we have analyzed in detail several aspects of the turduckening technique for 
evolving black holes. 

First we presented a detailed analytical study of the constraints propagation for a rather general 
family of BSSN-type formulations of the Einstein equations. We could appropriately identify a sub- 
family for which the constraints propagate within the light cone and give a rigorous justification of 
the turduckening procedure. At the same time, we showed that in other subfamilies the constraint 
violations do move superluminally. As a consequence, in those cases, smoothing the interior of 
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Figure 12: Difference between the I = 2, m = mode of r^4 for the turducken runs with n = 6 and n = 2 
(both with rt = 0.1) extracted at R = 30M (left plot). The right plots shows the same for the puncture run 
and the turducken with rt = 0.1 and n = 2. 



black holes will result in constraint violations that propagate to the outside. 

Through high-resolution spherically symmetric numerical simulations we analyzed in detail the 
behavior of the constraints and gauge modes. In particular, we confirmed to very high accuracy 
the predictions of our analytical study. We found that the numerical constraints are preserved 
outside the black hole when the theory predicts so, and are violated otherwise. We also found 
cases in which gauge modes do propagate superluminally, escaping from the black hole, and cases 
in which they are trapped in the inside. These differences in gauge modes have consequences for 
wave extraction, as discussed below. We also observed that, when the constraints are guaranteed 
to propagate within the light cone, the region of constraint violations inside the black hole shrinks 
with time, and that the same final stationary configuration seems to be approached, regardless of 
the details of the turduckening procedure. We also provided explanations for these features. 

Finally, we presented detailed three-dimensional simulations of single distorted black holes, 
comparing turduckened and puncture evolutions. We studied the effect that these different methods 
have on the coordinate conditions, constraint violations, and extracted gravitational waves. We 
found the waves to agree up to small but non-vanishing differences. Our convergence tests showed 
that those differences are not numerical artifacts but true features of the solution, caused by 
superluminal gauge modes escaping from the black hole. We also found that these differences in 
waveforms decay with increasing extraction radius. 
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